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Preface 


The  purpose  of  this  study  was  to  analyze  how  a  re¬ 
flected  wave  is  disturbed  as  it  propagates  from  a  reflecting 
surface  to  a  receiver  plane  via  a  thermally  bloomed  medium. 
Previous  studies  on  this  subject  have  been  done.  However, 
this  study  incorporates  more  sophisticated  models  for  the 
thermal  blooming  and  ray  propagation  portions  of  the 
analysis. 

Along  with  the  return  wave  analysis,  this  report 
includes  a  computer  program  developed  to  perform  a 
geometrical  ray  trace  through  the  thermally  bloomed  medium. 

The  program  is  thoroughly  documented  and  can  be  easily 
adapted  to  any  blooming  model.  In  its  present  form,  the  [  Acccr.?-' - 


program  is  adapted  for  use  with  the  Air  Force  Weaponsj 
Laboratory  blooming  model  known  as  PPiOPMD.  | 
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Abstract 

A  computer  code  was  developed  to  perform  a  geometrical 
ray  trace  of  a  light  ray  from  a  target  plane  to  a  receiver 
plane  via  a  thermally  bloomed  medium.  Thermal  blooming  was 
modeled  using  an  Air  Force  Weapons  Laboratory  computer  code 
called  PROPMD. 

Ray  traces  were  done  for  various  degrees  of  thermal 
blooming,  and  Strehl  ratio  calculations  were  made  for  each 
situation.  Also,  contour  plots  of  the  ray  optical  path 
deviations  across  the  receiver  plane  were  made. 

Results  showed  that  a  ray's  direction  of  propagation 
was  virtually  unaffected  by  refraction.  Disturbance  of  the 
ray  was  due  totally  to  the  effect  of  thermal  blooming  on  the 
optical  path  length  of  the  ray.  In  severe  blooming 
situations,  calculated  Strehl  ratios  were  less  than  0.1. 
However,  limiting  the  size  of  the  receiver  optics  to 
something  less  than  the  size  of  a  high-energy  beam  grid  led 
to  improved  Strehl  ratios. 
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RETURN  WAVE  ANALYSIS  IN  A  THERMALLY  BLOOMED  MFDIUM 


I  INTRODUCTION 


Background 

Propagation  of  a  high  power  laser  beam  through  an 
absorbing  medium  has  been  shown  to  be  a  nonlinear  process. 
As  the  beam  propagates,  some  of  its  energy  is  absorbed  by 
the  medium,  and  this  absorption  induces  temperature  changes 
within  the  medium.  The  temperature  changes  affect  the 
density  of  the  medium  which  in  turn  affect  the  refractive 
index  n  through  the  relationship 

n,  *  [«.- 1]  (i) 

where  nQ  is  the  refractive  index  of  the  quiescent  medium,  f o 
is  its  quiescent  gas  density,  n^  is  the  index  change,  andf. 
is  the  density  change.  The  relationship  between  temperature 
and  density  is  such  that  as  temperature  increases,  density 
(and  refractive  index)  decreases.  If  we  consider  the 
irradiance  profile  of  a  Gaussian  beam  where  the  center  of 
the  beam  is  the  hottest  area,  the  medium  resembles  a 
negative  lens  (see  Figure  1)  and  the  beam  is  defocused  (Ref 
1:27-19) • 


Bear*)  t>i  a  me  fer* 


Figure  1.  Heated  Medium  Resembles  Negative  Lens 


When  a  crosswind  (wind  traveling  in  a  direction 
perpendicular  to  that  of  the  laser  beam)  is  present,  heating 
across  the  beam  is  redistributed  and  the  refractive  index 
profile  for  the  beam  changes  (see  Figure  2) .  This  change  is 
due  to  the  cooling  effect  that  the  wind  has  on  the  beam.  As 
parcels  of  air  within  the  beam  are  heated  up,  they  are  sv.’ept 
out  of  the  beam  by  the  wind.  Since  the  downwind  portion  of 
the  beam  sees  less  of  a  cooling  effect  than  the  upwind 
portion,  the  index  profile  resembles  that  shown  in  Figure  2. 
In  addition  to  causing  the  beam  to  diverge  slightly,  this 
profile  also  causes  the  beam  to  bend  into  the  wind. 


WIMD - > 

Figure  2.  Refractive  Index  Profile  with  Cross wind 

From  the  standpoint  of  classical  geometrical  optics,  if 
the  refractive  index  of  any  point  within  the  beam  can  be 
determined,  then  a  Snell's  Law  analysis  of  a  light  ray 
traveling  inside  the  beam  can  be  performed.  Furthermore,  a 
calculation  of  the  optical  path  length  (OPL)  can  be  made  for 
any  particular  ray,  the  OPL  being  the  integral  of  the 
refractive  index  along  the  ray  path. 


With  the  help  of  an  Air  Force  Weapons  Laboratory  (AFWL) 
computer  code  called  PROPMD  and  the  basic  assumptions  listed 
later  in  this  section,  a  geometrical  optics  analysis  is 
possible. 

An  analysis  of  this  nature  is  of  interest  to  AFWL  s ^ nee 
low  power  lasers  (serving  as  tracking  lasers)  share  the  same 
aperture  as  and  travel  colinearly  with  high  power  laser 
beams.  If  it  were  possible  for  the  laser  operator  to  form 
an  image  of  the  target  by  using  the  return  wave  reflections 
from  the  target,  then  the  laser  operator  could  determine  not 
only  if  he  was  hitting  the  target  but  also  what  specific 
area  of  the  target  he  was  hitting. 

Problem 

The  problem  investigated  in  this  study  is  the  analysis 
in  the  receiver  plane  of  a  wave  that  has  emanated  from  a 
point  in  the  target  plane,  and  traveled  through  a  thermally 
bloomed  medium  to  the  receiver. 

The  analysis  does  not  include  the  effects  of 
diffraction  on  the  propagation.  Furthermore,  the  analysis 
is  related  only  to  the  return  wave  prior  to  its  entering  the 
receiver  imaging  optics. 
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Assumptions 

It  must  be  emphasized  that  the  accuracy  of  all  results 
of  this  research  are  dependent  upon  the  overall  accuracy  of 
the  PROPMD  computer  code  calculations.  Based  on  an  analysis 
of  PROPMD  calculations  and  on  certain  properties  of  the 
PROPMD  code,  the  following  assumptions  are  made: 

a.  The  value  of  the  index  of  refraction  for  any 
particular  x-y  coordinate  point  in  an  index  screen  (an  array 
of  index  values  in  a  plane  perpendicular  to  the  direction  of 
ray  propagation)  remains  constant  for  that  x-y  position  from 
a  z  coordinate  point  midway  between  that  index  screen  and 
the  preceding  index  screen  to  a  z  coordinate  point  midway 
between  that  index  screen  and  the  next  successive  index 
screen . 

b.  Based  on  the  first  assumption,  refraction  takes 
place  only  at  points  in  the  plane  midway  between  any  two 
index  screens. 

c.  The  index  of  refraction  value  used  in  calculating 
the  OPL  between  any  two  successive  index  screens  is  (through 
linear  interpolation)  the  average  value  of  the  two  indices. 

The  basic  assumptions  of  the  PROMPD  code  are: 

a.  The  beam  has  been  on  long  enough  to  establish  a 
steady-state  irradiance  distribution. 

b.  The  heat  transfer  process  within  the  beam  is 
instantaneous,  thus  producing  instantaneous  temperature 
changes  that  translate  to  instantaneous  refractive  index 
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changes . 


c.  The  high-energy  beam  has  an  unperturbed  irradiance 
profile  that  is  Gaussian  in  shape, 

d.  The  medium  through  which  the  beam  passes  j  s  free  of 
turbulence  and  is  flowing  at  a  constant  velocity  that  is 
much  less  than  the  speed  of  sound. 

The  wavelength  of  the  high-energy  beam  used  in  this 
study  is  10.6  micrometers. 

No  consideration  is  given  to  rays  that  reach  the  edge 
of  the  beam  grid  prior  to  reaching  the  receiver  plane.  The 
rays  are  allowed  to  leave  the  grid  and  are  not  j. ncluded  in 
calculations  made  in  the  receiver  plane. 

Development 

The  design  of  the  ray  trace  algorithm  is  presented  in 
chapter  II.  Special  attention  is  given  to  certain  portions 
of  the  computer  code  that  are  special  adaptations  of 
otherwise  general  knowledge. 

In  Chapter  III,  the  methods  of  anaylsis  are  discussed. 
The  various  results  of  the  analysis  are  compared  in  Chapter 
IV,  and  the  conclusions  and  recommendations  are  presented  in 
Chapter  V. 


II  COMPUTER  CODE  DEVELOPMENT 

The  most  difficult  and  time  consuming  portion  of  this 
study  has  been  the  development  and  validation  of  the 
computer  code  that  performs  the  geometrical  ray  trace. 
However,  before  the  ray  trace  program  could  be  used,  certain 
other  data  were  needed. 

PROPMD  Data 

E a  •  1  y  in  this  research  effort,  PROPMD  was  studied  in 
order  to  understand  the  mechanics  of  the  program  and  where 
certain  information  is  calculated  within  the  program.  A 
discussion  of  how  PROPMD  functions,  the  changes  made  to  it, 
and  a  sample  listing  of  the  input  parameters  used  can  be 
found  in  Appendix  B. 

Design  Intent 

The  ray  trace  program  called  TRACE  is  designed  to 
propagate  a  ray  from  a  specified  position  in  the  target 
plane  to  the  receiver  plane  through  a  medium  whose 
refractive  index  varies  as  a  function  of  position.  As  the 
ray  propagates,  its  initial  direction  and  velocity  are 
affected  by  the  variations  in  density  and  refractive  index 
from  point  to  point.  This  requires  a  three-dimensional  (3- 
D)  Snell's  Law  refraction  (see  Appendix  A  for  detailed 
development)  of  the  ray  at  points  along  the  propagation  path 
where  refractive  index  changes  are  significant. 
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Additionally,  a  cumulative  calculation  of  the  OPL  is  needed 
with 

pISc 

OPL  —  Zj  /  n6*)dr 

ISt-, 

where 

n  =  refractive  index 

r  =  the  vector  between  the  two  index  screens 
IS^  =  ith  index  screen 

As  stated  in  the  previous  chapter,  the  linear  interpolation 
between  planes  gives 

'AM -  [nui  +  vuJ/2  “  *U  (3) 

Thus,  the  form  of  the  OPL  equation 

rISL 

OPL  =  II J  rtiJr  =  2  »U  n  (4) 

ISc-, 

is  used  for  all  calculations  within  TRACE. 

Index  Screen  Locations 

TRACE  is  designed  to  be  given  the  location  of  the  index 


screens  along  the  propagation  path.  As  each  index  screen 
location  is  read,  the  appropriate  midplane  locations  are 
determined  for  later  calculations  of  refraction  and  OPL. 


Index  Interpolation 

As  a  ray  propagates  to  a  new  mid  plane,  its  angle  of 
incidence  (with  respect  to  the  surface  normal  of  the 
midplane)  is  calculated  as  well  as  its  new  x-y  coordinates. 
The  next  step  is  then  to  determine  the  value  of  n  at  that  x  - 
y  coordinate  in  the  index  plane.  This  determination  of  the 
value  is  based  on  a  linear  interpolation  subroutine  called 
ANDEX.  ANDEX  reads  in  a  square  matrix  representing  the 
index  distribution  across  the  index  plane.  Based  on  the 
spacing  between  array  values  (grid  spacing  value  determined 
by  PROPMD) ,  a  linear  interpolation  is  performed.  Since  the 
coordinate  system  of  PROP  HD  is  opposite  to  that  of  the  ray 
trace,  the  interpolation  begins  in  the  upper  right-hand 
corner  of  the  index  grid  where  n(1,1)  is  located. 

Optical  Path  Difference  ( OPD)  In  Receiver  Plane 

In  order  to  determine  what  changes  the  image  : n  the 
receiver  plane  has  undergone,  it  is  necessary  to  determine 
how  much  the  phase  front  of  the  actual  image  wave  differs 
•  from  that  of  a  perfectly  spherical  wave  that  has  propagated 
the  same  range  distance  through  an  unperturbed  medium.  This 
ideal  value  of  OPL  is  subtracted  from  the  cumulative  OPL 
value  for  each  ray.  To  find  the  true  CPD  within  the 
receiver  plane  for  each  ray,  ^  (see  Figure  3)  is  found 
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u  s  i  n  g 


(s) 


l  z  'l  X1  +  Yl  ♦  Zl  -  Z 

The  final  OPD  is  found  using 

OPD  ~  OPL  ~  S  rif  (&) 

where  n  ^  is  the  final  refractive  index  for  any  particular 
ray  (Ref  3:197).  Figure  3  assumes  that  nQ  is  the  final 
index  value  and  that  £  is  small  compared  to  z. 


Figure  3.  Geometry  of  OPD  Calculation. 

Strehl  Ratio  (Refs  4:151 3 ; 5 : 460  —  46 4 ) 

If  the  variations  in  OPD  across  the  receiver  plane  are 
small  enough  (on  the  order  of  0.5  wavelengths  or  less),  then 
the  Strehl  ratio  (a  ratio  of  the  peak  intensity  in  the  focal 
plane  of  a  disturbed  wave  front  to  the  peak  intensity  in  the 
focal  plane  of  an  ideal  undisturbed  wave  front)  can  be 
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approximated  by 


S+rekl  =  |  ~  £  OPO*  "  J  ^ 

where 

=  the  average  of  the  square  of  OPD 

__a 

QpD  =  the  square  of  the  average  OPD 
)y  =  the  wavelength 

Receiver  Plane  Mask j  ng 

In  order  to  analyze  how  the  Strehl  ratio  varies  as  the 
size  of  the  receiver  plane  being  analyzed  is  varied,  a 
parameter  call  MASK  is  input  to  TRACE.  This  parameter 
determines  what  radius  with  respect  to  the  center  of  the 
receiver  plane  is  analyzed  in  the  Strehl  ratio  calculations. 

Program  Verification 

The  validity  of  the  refraction  and  propagation  portions 
of  TRACE  were  verified  by  simulating  a  propagation  through 
three  phase  screens.  The  index  values  at  each  plane  were 
constants.  The  output  values  of  TRACE  were  checked  against 
hand  calculations  and  found  to  be  consistent. 

The  interpolation  portion  of  the  program  was  verified 
by  inputting  a  4x4  index  array  to  TRACE.  The  previously 
determined  indices  for  selected  points  on  the  index  grid 
were  then  compared  to  the  interpolated  values  from  TRACE. 
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Ill  ANALYSIS 


Since  there  are  so  many  parameters  in  the  thermal 
blooming  model  that  can  be  varied,  all  parameters  except 
power,  number  of  screens,  and  crosswind  velocity  are  kept 
constant  (see  Appendix  B  for  a  complete  listing  of  PROPMD 
parameters).  All  analysis  is  for  a  range  of  1000  meters, 
and  the  power  levels  studied  are  100,  10,000,  25,000, 
50,000,  75,000,  and  100,000  watts.  Analysis  of  these  power 
levels  is  done  for  crosswind  velocities  of  2.57  and  144. 
meters  per  second  for  both  5  and  10  screens. 

Due  to  the  problems  with  computer  access  and  job  turn 
around,  the  size  of  the  index  array  at  each  index  screen  is 
limited  to  64  x  64.  This  array  size  directly  affects  the 
distance  between  array  points,  and,  ultimately,  the 
resolution  of  index  changes. 

Receiver  Plane  Contours 

As  discussed  earlier,  a  high-energy  laser  beam  will 
bend  into  the  wind  when  a  cross wind  is  present,  and  the 
downwind  portion  of  the  beam  will  have  a  lower  refractive 
index  than  the  upwind  portion.  One  will  then  see  a  symmetry 
transverse  to  the  crosswind  direction  when  contour  plots  of 
the  OPD  in  the  receiver  plane  are  plotted  for  waves 
originating  from  an  on-axis  point  in  the  target  plane. 
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Strehl  Ratio 

The  Strehl  ratio  is  the  standard  of  measurement  for 
return  wave  quality  in  this  study.  Since  the  Strehl  ratio  is 
directly  related  to  phases  and  phases  are  directly  related 
to  OPD's,  the  Strehl  ratio  calculation  can  be  made  using  the 
final  OPD  values.  Furthermore,  it  might  be  possible  to 
determine  if  the  return  wave  over  a  certain  area  of  the 
receiver  plane  is  less  affected  than  that  of  another  area. 
Part  of  the  analysis  deals  with  a  masking  of  the  receiver 
plane  such  that  the  Strehl  calculations  include  only  these 
values  of  OPD  falling  within  a  previously  specified  radius 
of  the  center  of  the  receiver  plane. 

Location  Of  Blooming.  Effects 

In  an  attempt  to  determine  at  what  point  along  the 
propagation  path  the  thermal  blooming  has  its  greatest 
effect  on  the  return  wave,  various  numbers  of  index  screens 
are  skipped  before  refraction  takes  place.  Up  to  and 
including  the  screen  skipped,  the  OPL  calculation  uses  n  Q 
so  as  to  simulate  propagation  through  an  ambient  medium. 
The  Strehl  ratios  with  the  various  screens  skipped  can  then 
be  compared  to  determine  which  screens  contribute  the  most 
to  OPD. 
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IV  RESULTS 

The  findings  being  discussed  in  this  section  are  the 
result  of  analyzing  data  generated  through  the  following 
process : 

1.  Executing  the  PROP  MI)  program  using  a  selected 
combination  of  input  parameters. 

2.  Executing  the  DATA  program  to  translate  the  PROPMD 
program  output  into  the  required  sequence  of  index  screens. 

3.  Executing  the  TRACE  program  to  calculate  the  final 
x-y  positions  and  OPD  values  for  the  input  values  of  ray 
origin  and  launch  angle  representing  a  bundle  of  rays 
directed  from  the  target  back  to  the  receiver  (TRACE  uses 
the  index  screen  information  to  do  the  refraction  and  OPD 
calculations) . 

4.  Drawing  contour  plots  and  making  Strehl  ratio 
calculations  using  the  OPD  values  for  those  rays  that  have 
propagated  to  the  receiver  plane. 

The  following  discussion  attempts  to  explain  the 
findings  in  each  major  area  of  analysis.  Particularly 
significant  results  are  highlighted,  and  comparisons  are 
made  of  data  differences  associated  with  the  variation  of 
PROPMD  input  parameters. 

Rav  Refraction  and  Tilt 

In  this  portion  of  the  analysis,  a  contour  plot  of  OPD 
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values  for  each  ray  with  respect  to  that  ray’s  final  x-y 
position  in  the  receiver  plane  was  made.  As  expected,  these 
plots  showed  the  OI’D  to  be  symmetric  with  respect  to  the 
crosswind  flow  direction.  As  can  be  seen  in  Figures  4,  5, 
and  6;  the  contours  are  similiar  with  only  the  distance 
between  the  OPD  contour  lines  chancing.  From  a  total 
variation  in  OPD  of  0.5  waves  (waves  =  OPD/ X )  at  25,000 
watts,  we  see  a  variation  of  2.0  waves  at  100,000  watts. 

When  comparing  each  ray's  final  x-y  position  to  the 
projected  position  for  that  same  ray  had  it  traveled  through 
a  homogeneous  medium,  we  find  that  there  is  no  deviation  of 
the  rays.  This  comparison  implies  that,  for  the  beam  powers 
and  ranges  studied  herein,  ray  refraction  effects  are 
negligible . 

However,  there  is  a  nearly  linear  variation  of  OPD 
versus  distance  in  the  direction  of  wind  flow.  A  linear 
variation  of  phase  is  equivalent  to  tilt  of  the  wavefront, 
and  phase  can  be  directly  related  to  OPD.  The  tilt  is 
related  to  the  bending  into  the  wind  of  the  high  energy  beam 
that  disturbed  the  medium.  Wavefronts  propagating  through 
this  disturbed  medium  also  tend  to  bend.  This  bending 
causes  returning  wavefronts  to  be  tilted  with  respect  to  the 
receiver  plane. 

The  refraction  technique  used  in  this  study  does  not 
correct  the  ray  direction  to  account  for  the  presence  of 
wavefront  tilt,  so  a  rough  calculation  was  done  to  estimate 
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the  wavefront  tilt  for  a  beam  power  of  100,000  watts.  From 
Figure  6,  the  change  in  C'PD  from  the  center  of  the  grid  to 
the  right-hand  edge  gives  an  estimated  tilt  of  approximately 
2  microradians.  For  a  range  of  1000  meters,  this  tilt 
causes  a  ray  translation  of  2  millimeters  in  the  receiver 
plane.  Compared  to  the  un tilted  ray  translations  for  this 
beam  power,  the  tilt  translation  is  small  and  can  be 
ignored.  However,  '-he  combination  of  higher  beam  powers  and 
longer  propagation  ranges  could  produce  a  significant  amount 
of  translation. 

Strehl  Calculations 

Since  the  variation  in  OPD  for  Figure  6  is  large,  the 
Strehl  ratio  for  100,000  watts  will  be  very  low.  Thus,  any 
imaging  system  looking  at  a  grid  area  the  same  size  as  that, 
of  the  high  energy  beam  will  be  severely  affected. 

The  effects  of  looking  at  smaller  areas  within  the 
receiver  plane  proved  to  be  very  important.  In  general, 
the  analysis  showed  that  by  looking  at  smaller  areas  within 
receiver  plane,  one  could  calculate  a  better  Strehl  ratio. 
Figure  7  is  a  comparison  of  Strehl  ratios  for  a  crosswind  of 
144  m/s.  The  Strehl  ratios  for  a  crosswind  of  2.57  m/s  were 
all  zero  and  did  not  require  plotting.  Even  for  powers  of 
75,000  to  100,000  watts,  the  Strehl  ratio  is  better  when  a 
receiver  radius  of  0.1  meters  is  used,  and  the  ratio  starts 
to  fall  off  rapidly  as  the  radius  is  increased.  This 
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observation  fits  well  when  you  consider  Figure  6.  Though 
there  is  a  significant  amount  of  OPD  variation  across  the 
receiver  plane,  the  variation  over  small  portions  of  the 
plane  are  themselves  small.  However,  it  should  be  under¬ 
stood  that  limiting  the  aperture  radius  reduces  the  input 
signal.  Thus,  this  reduced  input  signal  can  lead  to  a 
reduced  signal-to- noise  ratio  and  no  significant  improvement 
in  the  quality  of  a  final  image. 

Location  Of  Bloom j.ng  Effects 

An  attempt  was  made  to  determine  which  scree  n  or 
screens  contribute  most  to  the  reduction  of  the  Strehl 
ratio.  Figures  8,  9,  and  10  are  plots  of  the  Strehl  ratio 
when  successive  index  screens  are  skipped  or  eliminated  from 
the  Strehl  ratio  calculation. 

The  plots  reveal  that  the  thermal  blooming  effects  for 
higher  beam  powers  are  concentrated  in  the  middle  index 
planes  (at  200  to  300  meters  from  the  target  plane).  Index 
planes  nearest  the  target  have  little  effect  on  the  Strehl 
ratio . 

This  observation  does  not  appear  to  agree  with 
published  findings.  It  has  been  experimentally  determined 
that  blooming  of  a  focused  beam  will  reach  its  maximum 
somewhere  near  the  focal  plane  (Ref  1:49).  The  value  of 
"somewhere"  is  not  strictly  defined,  so  blooming  in  a  plane 
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250  meters  from  the  target  might  be  feasible.  However,  one 
might  expect  the  most  severe  blooming  to  occur  closer  to  the 
focal  (target)  plane  for  a  range  of  1000  meters. 

One  explanation  for  the  apparent  disagreement  in 
blooming  locations  is  the  fact  that  the  launch  angles  of  the 
rays  from  the  target  plane  are  small  enough  to  keep  the  rays 
close  together  as  they  propagate  through  the  first  few  index 
screens.  Since  all  rays  are  launched  from  the  same  point  in 
the  target  plane,  they  all  see  essentially  the  same  index 
variation  through  the  first  few  index  screens.  For  this 
reason,  the  point  at  which  thermal  blooming  effects  begin  to 
become  noticeable  is  shifted  away  from  the  target  plane. 

Even  if  the  rays  are  not  kept  close  together,  we  would 
still  expect  to  see  the  same  type  curves  as  those  in  Figure 
8,  9,  and  10.  Once  the  high  energy  beam  has  been  thermally 
bloomed,  one  might  expect  to  see  less  significant  charges  in 
refractive  index  between  that  point  and  the  target  plane. 
Thus,  those  planes  closest  to  the  target  plane  would  have 
little  or  no  effect  on  the  Strehl  ratios. 

However,  it  must  be  pointed  out  that  the  accuracy  of 
this  particular  analysis  is  dependent  upon  the  accuracy  of 
the  information  determined  in  PR0PMD.  In  certain 
situations,  the  spacing  between  index  screens  could  greatly 
affect  the  resolution  of  index  variations  between  these 
index  screens.  For  beam  powers  greater  than  100,000  watts, 
we  would  expect  the  point  of  maximum  thermal  blooming  to 
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Figure  10.  Strehl  Ratio  vs  Screens  Skipped  (P=100,000  W) 
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shift  farther  away  from  the  target  plane.  However,  the 
PROPMD  program  spaces  the  planes  according  to  the  number  of 
index  screens  requested.  Thus,  for  higher  beam  powers,  the 
maximum  blooming  could  occur  in  areas  where  screen  spacing 
is  large. 
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V  CONCLUSIONS  AND  RECOMMENDATIONS 

The  objectives  of  this  study  were  to: 

a.  Develop  a  computer  code  based  on  geometrical 
optics  theory  for  propagating  rays  from  a  target  plane  to  a 
receiver  plane  via  a  thermally  bloomed  medium. 

b.  Analyze  the  effects  on  the  rays  of  propagation 
through  a  thermally  bloomed  medium. 

Both  of  these  objectives  were  met,  and  a  listing  of  the 
developed  computer  code  was  included  in  Appendix  C. 

The  analysis  showed  that  for  the  combination  of  high 
beam  powers  and  low  crosswind  speeds,  the  return  wave  in  the 
plane  of  the  receiver  was  degraded  to  the  point  that  the 
Strehl  ratio  was  always  less  than  0.1.  Additionally, 
limiting  the  size  of  the  receiver  aperture  did  not  improve 
the  Strehl  ratio  no  matter  how  small  the  aperture  was. 

For  a  crosswir.d  of  1 H  *1  m  /  s ,  Strehl  ratios  could  be 
improved  by  using  smaller  receiver  apertures.  However,  one 
should  be  aware  that  tradeoffs  must  be  made  with  respect  to 
the  resulting  input  signal  reduction  and  possible  reduction 
of  the  signal-to-noise  ratio. 

One  significant  result  of  this  study  was  that  ray 
refraction  caused  by  thermal  blooming  was  negligible.  Thus, 
the  rays  were  disturbed  only  with  respect  to  their  OPD 
values.  However,  there  was  a  tilt  in  the  returning 
wavefront.  Since  the  propagating  high  energy  beam  was  bent 
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into  the  wind,  wavefronts  traveling  the  same  path  were  also 
bent.  This  bending  caused  returning  wavefronts  to  be  tilted 
with  respect  to  the  receiver  plane. 

Because  light  rays  travel  perpendicular  to  a  wavefront, 
wavefront  tilt  implies  that  the  light  rays  are  bent.  In 
this  study,  the  affect  of  wavefront  tilt  on  ray  refraction 
was  ignored,  and  an  estimate  on  the  amount  of  tilt  for  the 
100,000  watt  case  was  not  considered  large  enough  to 
invalidate  the  results  of  this  study. 

However,  the  effects  of  more  significant  amounts  of 
tilt  must  be  considered  since  the  TRACE  program  does  not 
accurately  model  the  refraction  of  propagating  rays.  If 
this  tilt  is  large  enough  in  comparison  with  the  original 
launch  angle  of  an  individual  ray,  then  it  can  also  have  a 
significant  effect  on  the  OPD  calculations.  Even  if  the 
tilt  is  not  significant  at  a  range  of  1000  meters  and  a  beam 
power  of  100,000  watts,  it  must  be  considered  for  different 
ranges  and  beam  powers.  Higher  beam  powers  will  cause  more 
beam  bending,  and  longer  ranges  will  amplify  the  translation 
effects  of  tilt  in  the  receiver  plane. 

Depending  on  the  accuracy  of  the  PR0PMD  calculations, 
the  spacing  of  the  index  grid  plane  can  have  a  significant 
effect  on  the  resolution  of  index  changes  between  grid 
points  within  an  index  plane.  Therefore,  in  order  to 
achieve  the  best  possible  resolution,  the  optimum  grid 
spacing  (optimum  index  array  size)  should  be  determined. 
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As  a  starting  point,  a  linear  interpoi ation  scheme  was 
used  to  determine  the  index  variations  between  two  index 
planes.  Even  though  the  process  of  thermal  blooming  is 
nonlinear  in  nature,  a  linear  interpolation  scheme  is 
appropriate  if  PROPMD  spaces  the  index  screens  properly. 
However,  the  results  of  determining  the  location  of  thermal 
blooming  effects  indicates  that  a  higher  order  interpolation 
scheme  would  be  more  accurate. 

Based  on  the  above  conclusions,  the  following 
recommendations  are  made: 

1.  The  imaging  optics  should  have  the  smallest  field- 
of-view  possible  for  the  given  minimum  system  requirements 
on  input  signal  strength  and  signal-to-noi se  ratio. 

2.  Further  study  should  be  done  at  higher  beam  powers 
and  longer  ranges  with  effects  of  wavefront  tilt  included. 

3.  Further  study  should  be  done  to  determine  the 
optimum  grid  spacing  necessary  for  accurate  resolution  of 
index  variations. 

4.  Investigation  into  the  development  of  a  higher 
order  model  for  index  variations  between  index  planes  should 
be  made. 
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APPENDIX  A 


Lunebur g  1  s  3-D  Re fr act j  on  (Ref  6: OR -6 6) 

The  Lunebur g  3-D  refraction  technique  is  generalized 
for  any  s h a p e cl  surface  of  interface  between  two  media  of 
different  refractive  indices  (see  Figure  11). 


Figure  11.  Geometrical  Representation 

M  i s  a  unit  vector  in  the  direction  of  the  surface  normal  of 
S.  T  represents  a  unit  vector  in  the  direction  of  the 
incident  ray,  and  T 1  is  the  unit  vector  of  the  refracted 

ray . 

Luneburg’s  analysis  develops  the  equation 

Yi’T*  -  nT  s  PM  (6) 

where  P  is  a  scale  factor. 

Since  it  it  known  that  the  refracted  ray  leaves  the 
surface  in  a  plane  formed  by  the  incident  ray  and  the 


surface  normal  H,  it  follows  that  the  following  vector 


equation  is  true  .  . 

n'  (t'/M)  =  n  (TxM) 


The  lengths  of  the  two  vectors  in  equation  9  yield  Snell's 
1  aw 

sm  4*  -  sm  4*  (l°) 

where 


<$  is  the  angle  of  incidence 
<J)f  is  the  angle  of  refraction 
The  factor  f7  is  found  by  forming  the  scale r  product 


r 

=  n'  (t'-m)  ~  n  ('T-M) 

(H) 

or 

r 

=  Yl'  COS  4 1  "  Fl  COS  (j‘ 

02) 

The 

equation  of  the  surface  normal  M  in  this  study  is 

M 

A  A  i  A 

rox  +  OY-+  1  2 

(13) 

Let 

T 

=  <,  S  +  Y,  y  ♦  *,  S' 

(H-) 

and 

T* 

A  A.  —  ^ 

=  ^  x  ♦  Y  +  Z 

(is? 

Then 

equations  8, 

1 3 j  and  15  combine  to  give 

Y\‘  (*2  X  +  Yz  ’Y  ♦  2zl  )-  r\  (x.x-r  Y,Y  t2yz) 


=  (ox  +  ov-*  (lb) 
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Equation  16  can  be  rearranged  to  yield 


=  (vW  X, 

0 n 

'fi  ~(%‘)  V| 

(18 

i%-  (r  +  nz.)/V 

(ii. 

Since  T,  T1,  and  M  are  unit  vectors, 

j'x  M  -  Sm  (j/ 

(20) 

T  x  ~  S  m  <j> 

(21) 

Therefore,  equation  9  gives 

<j>*  c  ARCSh/  [(*/,)*" 

♦] 

(22) 

Combining  equations  12  and  19  yields 

Zz  -  (p  +  ii  z, ) /V  = 

COS  <J>* 

(23) 

and  the  refracted  unit  vector  T'  is  found. 
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APPENDIX  B 

A  Thermal  Blooming  Model 
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APPENDIX  B 


PROPMD  A  Thermal  Blooming  Mode] 

PROPMD  is  a  very  sophisticated  computer  code  written  by 
Dr.  Charles  B.  H  o  g  g  e  of  the  Air  Force  Weapons  Laboratory. 
The  program  has  been  designed  to  allow  great  latitude  in 
modeling  the  high-energy  beam  as  well  as  the  effects  of 
atmospheric  absorption  and  turbulence  on  the  beam’s 
propagation.  PROPMD  can  simulate  the  firing  of  a  laser  beam 
from  either  a  stationary  or  moving  platform.  A  crosswind  at 
any  speed  and  at  any  angle  with  respect  to  the  beam 
propagation  path  can  be  selected,  and  the  beam  itself  can 
travel  at  any  angle  with  respect  to  horizontal  (i.e.  an 
airborne  platform  with  the  beam  aimed  downward  and  to  the 
rear  of  the  aircraft).  There  are  input  parameters  that  can 
simulate  a  particular  atmospheric  turbulence  situation.  The 
beam  characteristics  (power,  wavelength,  beam  size,  pulse 
length,  etc.)  are  input  along  with  appropriate  values  of  the 
absorption  coefficients  and  relaxation  times  for  carbon 
dioxide  (C02)  and  water  in  order  to  model  the  absorption  of 
the  beam,  particularly  the  C02  laser  beam  of  wavelength  10.6 
micrometers . 

The  actual  location  of  phase  (index)  screens  can  either 
be  input  by  the  user  or  left  to  the  program  to  calculate 
internally  depending  on  how  many  screens  are  requested.  For 
5  arid  10  screens,  the  internally  generated  s pacings  are  as 
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shown  in  Figure  12.  Since  the  thermal  blooming  has  its 
greatest  effect  near  the  target  plane,  the  screens  are 
necessarily  close  together  in  that  area  of  the  propagation 
path . 

One  important  feature  internal  to  the  program  is  its 
compression  of  the  x-y  grid  as  the  beam  propagates  toward 
the  target.  The  overall  grid  size  of  the  propagating  beam 
is  reduced  to  simulate  focusing  of  the  beam.  Consequently, 
the  spacing  between  index  array  elements  is  reduced. 
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Figure  12.  Phase  Screen  Locations  In  FRGPMD 
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The  determination  of  thermal  blooming  effects  at  each 
screen  is  actually  a  calculation  of  phase  perturbations 
across  the  beam  grid.  These  perturbations  are  directly 
related  to  the  refractive  index  n  and  can  be  easily 
converted . 

For  persons  who  wish  to  use  PROPMD  in  conjunction  with 
TRACE,  the  following  additions  to  PROPMD  are  needed: 

1.  In  the  program  card  of  DEVICE,  TAPE22  should  be 
added . 

2.  After  line  131  of  DEVICE,  add: 

WRITE(22)  CITOT,  WAVE,  RANGE,  PETAM,  NSTP(I) 

3.  After  line  97  of  PROPMD,  add: 

DXPX  =  PX(2)  -  PX(1) 

WRITEC22)  DXPX,  NLM 

4.  After  line  995  of  PROPMD,  add: 

WRITEC22)  K,  AZ 

5.  After  line  49  of  CVTHM,  add: 

DXPX  =  PX ( 2 )  -  PXC1) 

WRITE(22)  DXPX 

DO  969  11=1,  NLM 

969  WRITEC22)  (UUS(I  I ,  J  J ) ,  J  J  =  1  ,  NLM ) 

After  a  successful  run  of  PROPMD,  the  data  written  on 
tape  22  must  still  be  manipulated.  The  data  is  oriented  to 
propagation  towards  the  target.  Whereas,  the  ray  trace  is 
away  from  the  target.  Thus,  the  order  of  the  screens  must 
be  reversed. 
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In  addition,  the  following  sample  operation  must  be 

performed  on  each  screen's  array  (see  DATA  in  Appendix  C): 

DO  100  1=1,  NLM 
A(  1 ,  1  )  =  1 .  +  PETAM 
DO  100  J  =  2,  NLM 

A  ( I ,  J  )  =  A(I,  J)/(2PI/WAVE)  +  (1.  +  PETAM) 

100  CONTINUE 

where  2PJ.  =  6.283  1  8. 

Using  the  values  of  RANGE  and  AZ,  the  spacing  between 
screens  can  be  calculated. 

The  following  is  a  list  of  sample  PROPMD  input 
parameters  used  in  this  study: 


CITOT 

= 

.1  E-6 

RAD  IN 

= 

0. 

WAVE 

= 

.106  E-4 

RADOUT 

= 

0.29  E  +  0 

VEL 

= 

.144  E  +  3 

BEAMSZ 

r 

.29  E+1 

RANGE 

= 

.1  E  +  4 

BETA 

r 

.2  E  +  1 

FL 

= 

.1  E  +  4 

TURBCST 

s 

.1  E-23 

PALPWAT 

= 

.19  E-4 

GRNDLEV 

= 

.54  E  +  4 

PALPC02 

= 

.627  E-4 

RATE 

r 

.1  E+1 

PTAU 

= 

.85  E-5 

TPULSE 

= 

.1  E  +  1 

PETAM 

= 

.272  E-3 

ADAP 

= 

.1  E  +  1 

PPO 

= 

.1013  E+4 

NRSS 

= 

1 

AZANGLE 

= 

-.9  E+2 

NSTPSIZ 

= 

0 

ELANGLE 

= 

0. 

GARE0 

- 

0 .  ,  0 . 

TMSDFLT 

= 

.1  E  +  1 

PSCAT 

= 

.1  E- 1  3 

RATIO 

= 

.1  E  +  1 

RHO 

= 

0. 

C00RUPD 

- 

.15  E  +  0 

BLOSS 

= 

.1  E+1 

THETAJ 

- 

0. 

SOASP 

- 

0. 
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SOOL 


0 


ALTTRNS 


SEG  =  . 1 E+1  ,  99*0.  NSTP 

NSEGMAX  =  1  WAVESBQ 

COREL  =  0 


APPENDIX  C 

Listings  for  TRACE  and.  DATA 


HI 


FRO'  -F  hM  1  Fhi'E  1  INPUT  >  OUTPUT  ,  THr'cV;.::H!  iTFi  IT  .  Th!:'f 

D  I  MEN*  I  On  FAY  •  1  0  0. 7  •  <  AULA  I  T  1  1  0  0  .  Hr  i  FrtP  i  n  || , 

D I  HE  f  C I  Cj f  <  ' !  1  64  >  «  YY  <  t-4  >  <  Hit t  DF ' >.  64  *  64 ) 

D I  HEM  1 1  ON  X  0  •  1  0  0 )  <  YD  •  1  0  0  >  .  A l. P : □  .  1  Ij  0 .  ,  H L P V  □  (  1  0  n ':- 

DIMENSION  GPD I FF  (  1  0 0> 


1  .  T H P E '? <  T 6 PE 7> 
MF LAN  1  00) 


■  ^  >  *  t>  ■*>  ** 


•»  THIS  PPD'7-Fhm  Ii:  DESIGNED  1  0  PE  USED  I H  C  DM  JUMP  T  I  DM  UITH  THE 

♦  5  PLI.-MI'  CODE  DEVELOPED  In'  6-I.IL-  nFTO.  FRGEND  CALCULATES  PHASE 

♦  iCPEEM.:  THAT  CAM  BE  t'E'i  I F  ULATED  TO  PRODUCE  A  SlFEEN  Dr 

♦  refractive  indexes  that  hFE  THEM  PE  PH  BY  THi:  PROGRAM. 

♦  THI  S-  PPD6P Fit'l  IS  BE.  IGNEP  TG  ChLl ULhTE  m  gECXCTF’ IlhL  F'Ay 

♦  TRACE  DF  THE  RETURN  PEFLECT  I  DM  FF  DM  A  TAPPET  ILLUMINATED 

♦  BY  BOTH  A  H I  AH  tNEF'GY  LA.  r.F'  BFhM  AND  H  LCU  F'0'"ic.F  TF'ACKER 

♦  BLAH  THAT  APE  I J  :■  I  f i u  A  ... H A P‘ E  Ii  A  R  L  F'  A  T IJ  F'  E  • 

♦  THIS  PROGRAM  1,1  Hi  DEVELOPED  AT  THE  A  IP  FDFCE  INSTITUTE  OF 

♦  TECHMOLDI3Y  By'  1LT  MICHAEL  T.  BAiEP.  IT  I, IAS  l,iP  I TTEH  IM 

♦  JULY  1  'A :•  0  AS  RAPT  DF  HIS  THESIS  RESEARCH. 

♦  4  >  4  * ♦*  4 ♦ 


♦  PAY  ( 1  .  1 ')  =  XO  (METERS') 

♦  RAY  •:  1  »  S' )  =  ALF'HA  X 

♦  PAY  C 1  *  3)  =  VO  (METERS') 

♦  F' A V  ( 1  i  4 )  —  ALPHA  Y 

♦  FAY < 1 i S'  =  ALPHA  Z 

♦  RAY  (.  1 1 6 )  =  OPTICAL  PATH  lF.NGTH  OF  THE 

START  I  MG  1,1 1 TH  OPB  =  0  FO 

♦  PRY  <  1  >  7')  =  REFRACTIVE  INDEX  (INTERPOLATED  AT  EACH 

♦  ALPHA  X«  Y.  AMP  Z  =  DIRECTION  CDS  I  MSS  CF  THE  RAY 

♦  ZPLAST  (ZPMEI.il )  AMD  ZSLH ST  ■  Z  SMFM",  =  OLD  (MEM')  DISTANCES  OF 


FAY  -CUMULATIVE 


START  I  hi? 


•DIMT  '•  METERS') 
FLAME 

THF 


MIDPlhMSS  AMD 

SCREENS  ALONG  Z  AXIS  (METERS') 

♦ 

MIC 

NUMBER 

OF  DELTA  Z  IM 

F  E M c  M T  S  1 1 1  F  F  OP  A6H  i  I  u M 

DELZ 

= 

Z  AXIS 

INCREMENT  BTI.l 

M  SCREENS  •  METERS') 

♦ 

TAPE'? 

THE  FILE  CONTAINING 

INFO  FOP  PLOTTING 

♦ 

TAPES'  1 

= 

THE  FILE  CONTAINS  IM 

F 0  FROM  A  FFORMD  RIJN. 

♦ 

MASK 

= 

p  l£  L»  ij  i.H  ; 

R‘ AB  T  Cl  _  OF  Pi.: 

v'R  PLA.ME  USED  FOR  STREHL 

CALC. 

10 

= 

STPflHL 

PhD I U I  CJ F  0.  1 

METERS  . 

£0 

= 

STF'EHL 

FA DIUS  0-  0.£ 

METERS. 

3  0 

= 

STREHL 

RADIUS  OF  fj.3 

METERS . 

4  0 

=r 

S T F  EHL 

RADIUS  OF  0.4 

METEPS. 

50 

= 

STREHL 

RADIUS  OF  0.5 

METERS  . 

55 

= 

STREHL 

F  A  I'  I  IJ  S  0  F  0 . 5 

5  METERS. 

♦♦*HDTE : 

M H  \  k  H P E H  f'l H !  I  M  I  3 cl1 

FOR  RANGE  OF  1000  METERS 

amalys i 

S  OF  OTHER  RANGES  MILL  F ECU  IRE  CHANGES 

IN 

THE  MAS 

LI  MG  LDiiIC-AL 

IF  STATEMENTS. 

♦ 

HFLA6 

MO.  OF 

SCREENS  TO  BE 

SKIPPED  IN  ANALYSIS 

U2 


♦  NF'PIMT 

0 

1 

c 

-1 

♦  N  F*  H  t 

♦  ML  M 

♦  NOTE: 


=  program  ppimt  cations 

=  WRITES  -FLEET E D  DATA  TO  TAPE  POP  LATEP  F'LDTT  IMG. 

•TAPES'  must  be  MADE  FEPMANENT  file. 

=  PRINTED  DUTPI.IT  duly. 

=  PR1MTED  OUTPUT  PLUS  WRITING  TO  TAPES'. 

=  ppimt:  f  P'DFtTi-  vhl.ijf:*  :tpehl  pat  id.  may  rev  i  at  i  dm 

=  NUMBER  DP  h  h  t  _>  TD  BE  PPCPAoriTED 

--  the  dimems  i cm  dp  imtey  affay:  fpdm  pfopmd. 
the  di men: ion:  dp  vv ,  amd  a  index  apfays 

MUST  BE  THE  _hmE  h.  THc.  ML  M  PAPFiMt.TEF'  FEnD 
IH  FPDM  TAP  Eel  FILE. 


♦  MOTE:  EFiCH  PA'.-’  IS  PEhIi  IM  DM  h  SEPARATE  CARD  IM 

SEQUENT  I  ML  OF'LEP  (F'h'i  •  1  *  1  >  *  .  .  .  .  .  PAY  *  1 * 7 :•  1  . 

♦  ♦♦MOTE:  PAY  >  3  MUST  HAVE  ALL  ELEMENTS  ZEPD  TD  MODEL  A  PAY 

STRAIGHT  DOWN  THE  MIDDLE  CF  THE  BEAM 

♦  ♦♦♦♦♦  ♦♦♦♦♦♦♦♦THI  --  P  A  i  fT i J T  BE  THE  F  I P' S' T  P A  r  I M  F’  !J  l  ^  *************  *• 

♦  A INDEX  =  THE  MATRIX  OF  H-S*Y>  READ  IN  FOR  EACH  SCREEN 

♦  A  S  D  P  D  =  A  V E  R AGE  OF  S  O U A  R E S  OF  G F D "  S  F f  P  OSS  OUTPUT  FLAME 

♦  SAOPD  =  SAUAPE  OF  AVERAGE  DP  CF'D'E  AC  FOSS  OUTPUT  PLANE 

♦  DELI  A  =  EXTRA  LENGTH  ERC  H  fhV  TRA'-'EL:  FP'OM  STRAIGHT  THPI.I  F  AY 

♦  DF'DIFF  -  DF'D  DIFFEFEMCE  FDR  EACH  FAY  AC  PC :  !•  THE  OUTPUT  PLANE 

♦  RMS U  =  VARIANCE  C-F  DF'D  DIFFERENCE  ACROSS  OUTPUT  PLANE 

♦  STREHL  =  sTFEHL  FAT  ID  FOR  THE  FAY  TRACE  —  VALID  WHEN  THE 

ALPHA  Z  IS  VERY  SM  ALL  — 


<►«  « ■ 


INITIALIZE  SOME  VARIABLES 


ZTOTAL=0. 
ZF'LhST = 0 . 
ZPNEW=0. 
ZSLAST-O. 
ZSNEW*0. 


♦  READ  IN  RAY  INFORMATION 

♦  ♦♦♦♦♦♦♦♦»♦  ♦^♦♦♦♦♦♦<.»*<>^ 


READ*  *  NR  AY*  NPPINT*  MASK*  NFL.  AG 
nflag=nelagm 
DO  10  1  =  1*  NR AY 
1  0  READ  7  *  *  )  '  PAY  <  I  *  J  >  »  J=  1  >  7  > 


«■  «•  +  (,  <r  «►  » 


♦  FIND  ALPHA  Z  FDR  INPUT  VALUES  DF  ALPHA  X  AND  ALPHA  Y 

I'D  11  I  ~  1 *  N R h Y 
'  .0  1  i  '  —  F  H  i  1  I  «  1  > 

I  L)  1  I  ■'  —  FA'i  <  1  *  j  1 

HLP  ■  I  1  —  PAY  -  1 , c ■ 

HI. PYl  I  1  =  F  AY  ■.  I  *  4  > 

C  hL  L  A;.  F  Z  •  PAY  '  I «  c  .  F  AY  >4 1 «  4  :•  *  PAY  ■:  1 . 5  :>  > 


■  fr  <*  ♦  + 


♦  READ  1 M  INFOFTlAT I DIjhl  VALUE’  DP  PRDPMH  CODE  THAT 

♦  CREATED  THE  INDEX  MATRICES. 

READ  •:  c  r>  PDMEP.  HAVE*  RANCH*  RET  AM 


NZEPD=1.  +  PETAM 
RANGE S-RANGE 


♦  overall  do  loop  fop  propagation  increments 


*♦*>  t>#  <><• 


'*■■>00 


AI  riEAL  =  PANGE«’NZERD 
READ  •'£  1'  NIC*  NLM 
NDZ 1 =  NDZ  +  1 
DO  1000  k'=  1  *  NDZ  1 
READ  •:  £  1  >  DELZ 
ZTOTAL=ZTOTAL  +  DELZ 


*  FIRST  TIME  THRU*  THE  INDEX  VALUES  IN  TARGET  PLANE  AP 
IF  <DELZ.  LE.  0.  0001:*  GO  TO  100 


♦  FIND  LOCATIONS  OF  INTERMEDIATE  SCREENS  AND  MIDPLANES 

CALL  ZP’R'OP  1  ZsLFi .  T  *  Z SNEIil *  DELZ*  ZFHEUs 


FOUND 


«►  «*  *  ♦  ♦  <•  <V  < 


♦  DD  LOOP  TD  PROPAGATE  ALL  PAYG  THRU  DELTA 


1  fm  CCMTINiJF 

r  j  p  * 1  *  i  i  —  i  %  f  ip  *  - 1 

hM^M .  I  •  I  • - K mV  *  I  %  »  * 

1 1'  •  i::i LZ.  L  t  .  0  .  0  0  m  1  •  GG  in  11  0 
1 F  »■  F  mV  <  1*7’*  .  LT  .  0 .  ••  Gu  1  □  5 0  0 


*  ^  V  <►  ' 


♦  First*  ANGLE  Lc  INCIDENCE 


C  M L  L  Hf  li.-L E  1  F  h'i  •  I « U  '•  »  THET Hi* 

♦  FI  HP  X  ::•  'V  C  DOP'D 1 1  IhTES  AT  HE  ill  M I  DPLAHE 


XL  A  ST  a  PAY'I.n 
YLAST  —  PAY  1  I  <  z*’> 

CALL  XYMEI.1 1  P’  A ' i'  i  I  »  1  )  «  PhV  (  1 1  c' *  *  PAY  •’  I «  o  >  <  P'AY'  1  I  » 
$  ZPNEl.il  <  PhY  >.  I  *  3  1  '  PhY  ■  I  *  1  1  *  F  h1/  <  I  *  3  *  > 


><>» 


♦  FIT  THIS  POINT  ME  MUST  INTERPOLATE  HEW  N  •: X 


110  CONTINUE 

IFCI.GT.  1*  GO  TD  ISO 
READ  '■  c'  1 >  NFL HUE  •  DELXY 
IiO  ISO  1I=1«NLN 

ISO  READ  'cl  )  < ft  INDEX  <  1 1  >  .J  J  >  %  *  NLMT 

130  GOUT  I  HUE 

CALL  AH  LEX  1  DELXY*  PAY  *  I  •  1>  <  RAY  <  I  ,  3')  »  HLM  »  RAY  >:  I « 
■l  A  INDEX*  XX.  YY«  I  > 

IF  •CRAY  <  I  «  T)  .  bT  •  0.  >  bO  TO  1  4 0 


INCREMENT 


’>  <  ZPLAST « 


Y> 


4»  <'  <>  <'  <*  ♦  <*  <* 


«-  IF  THE  IHDEO  IS  FOUND  TO  I.E  NEGATIVE*  THE  RAY  IS  OFF  THE 
«►  GR 1 1'  Mil  D  HO  HOPE  CALCULA1  ION.  HP  E  HEEDED  FOP  TH I  F'A  r . 

HPLAN'I'=  HP-LAME 
GO  TO  O'O 
140  CONTINUE 

♦  IF  NFLAG  IS  SET*  THEN  HFLHG  -  1  SCREENS  MILL  BE  SLIPPED  IM 

«>  THE  F  mV  TRACE  BEFORE  REF  FACT  I  ON  AMD  OF  D  CALCS  hFE  MADE 

1 F  *1 NF*  LAG  •  gE  •  t  )  HNL  mi  _  T  1  I  *  =  PH  i’  i  I '  §  ) 

IF  t'DELZ.  LE.  0.  0001  >  GO  TO  500 

♦  FIHD  RVEFHGE  I  HD EX  BETMEEH  H I BPLHHES 

CALL  AVGH  HlNLAsT  '■  I )  *  RAY  <I»?)  *  ANBAR  c  I )  > 

♦  <>  O  <-  <>  O  <v  «-«-*<►  ^  <v  <>  <.  O  C-  «.«.«►<>  <r<» 

♦  Finn  THE  OPD  BETHEEH  MI  BP  LANES 

IF  1K.LE.  NFLAG')  ANBAR  •:  I  :•  =NZERO 

CALL  DPL  '-PA  v  <  I »  & )  «  KLA-i-T «  PAY  •:  I »  1 )  »  YLh  .  T  *  PAY  ( I  j  3 «  ZF'LAST  «  2PNEW« 
i  ANBAR  ■:'!>> 


♦  BRING  THE  OPD  CALCULATION  UP-TO-DATE 

X=F'AY  '•  I *  1 ) 

Y=RAY  I ,  S') 


«•  r> 


♦  FIND  REFRACTION  ANGLE  AND  HEM  DIRECTION  COSINES 


•.-ML L  P'P  UN I  TV  (RAY  1  I  •  cl  ■  '  r-nv  1  1 « *1  >  «  F‘i lY  *.  1 «  5  >  »  At U.  A  I- 7  ••  I  ■'  «  PAY  1  1  «  7  >  • 
IT HE! Ml  ■ 

S o m  r UN t i •  me 

IP  •  >  MD71  :•  GO  TG  1  000 

Z;:  L !  .  Z Pf«E l.l 

;• ;  mem 

1  00  0  i' .  ;  r 


*  (r  «>$•<>  ^  <  t <►  i-  <•  O  O  <•<•<»<>«•<.  t'  <■  i  «►  «. 


♦  T Oil  E  CARE  DP  LAST  PROPAGATION  FROM  MIDPLANE  TO  TARGET  PL 


ZPL.Pi*  T  =  ZPMEW 
ZPNEM  =  Z'MEM 

DG  C  l.l  l.l  l.l  I  =  1  '  t  1  H  r!  '1 

IF  cP’PiY  i"I  t  ?<  .  LT.  0.  GO  TO  £000 
ANLAIT  I  ">  -PhY  d  <  7  > 

L  PlLL  PiNmLE  '  F'AY  •'  I  5  0  •'  ’  THE  TP)  I  > 

L H  .  1  —  PhV  (  ]  «  1  .' 

YLAS  I  ~  Phi' '  1 1 

CALL  IKYNEM  '•  PhY  *:  I «  1  :•  *  FPiY‘  <  I  <  £7  <  P  hY  <:  I «  3 *  P  hY  1  I  •  -?  ’>  <  ZF-'LH  I T 
‘I- *  iPMnlil «  F'AY'  >  I  i  G  1  '  PhY  '•  I  *  1 7  *  P'h'i  '  ]  <  .7  i  7 
CALL  OPL  1  RAY  t  I »  G  >  « ILL  Pi  ;T  <  Pm  P  <  1 «  1  j  »  Y'LAST  >  RAY  i  I  »  .S  >  <  ZF’LAS  T  «  ZPNG 
P'  A  V  1  I  .  7  ■'  •' 

£  000  CONT I  Mi.iE 

DO  £  100  I  —  1  >  t  iPPiY 
PAY  <  I  »  £7  =PhY  *:  I  *  G 7  -  Pi  I  DEAL 
£100  CONTI  HUE 

WAVEP=WAVE 
WAVE=  WAVE--'  1 .  OE-G 


♦  FIND  STPEHL  PATIO  AND  MAX  DEVIATION  ACROSS  OUTPUT  PLANE 


NSUB= 0 

DO  3000  1  =  1'  tiPAY 

DELT A  =  ii.;PT  1  P  A  v’  *  I  .  1 1  ♦PAY-1  1  I  *  1  .*  -fP’Pi'r'  1  I  *  ^P'Pi’v'  *'  1  *  Y 1  +PANhE.  ♦•'-’AMcE 

■I.ES 

IF  CRAY  ( I  .  ?':•  .  LT  .  0.  •'  N *  IJD=MS UB+  1 
OPD I FF  ■:  I  ■'  —  P  PiY  t 1 «  F.  <  — P  AY  *1  ?  7  *  ^DELl  A 
3000  IF  'PAY  .'I.  7  >  .LT.  0.7  OF  D  IFF  1 1  7  =  0 . 

NDI V=MPAY  -  NIUE 
SUMS 0=0. 

SUM=0. 


‘1 


L 


L=  0 

I'D  4  0  0  0  I 

•  1 <  HP AY 

RR---5-ART  '  f 

.  i  .  *  I  .  1  ..  ♦I-’A  Y  *. 

1  .  1  1  +  P  A  Y  I  .  3 

<*  PAY  •:  I 

I F  • :  •  i . r: '  i  . 

EC. 10'. AND. ' 

F  P  .  'j.  T  .  .  1  U  1  1 

GO 

TO 

IF  '  ■  r  •'< ;  •  . 

F  A .  3  0  >  .  A  r  <  r> .  >. 

P  P .  G  T  .  .  cl.'  i  •  1  )  :■ 

GO 

TO 

IF  .  ■  "XI  . 

r  i.1  •  *  i  1  •  h  *  i  ,  • 

<•  P'  .  G  T  .  .  U  1  1  .* 

GO 

TO 

IFF-  ■ ; .  . 

EA. 4  0'.  AND.  • 

R  *•  .  T  .  .  4  i.t  1  ■  > 

GO 

TO 

IF  (  *. !'.  Xi. 

FA.  50  '  .  AND.  . 

P'P  .  GT  .  .  5  C'  1  :> 

GO 

TO 

I F  •  •-  A  A ;  i  . 

L=L+I 

EC1.  55  1  .AND.  1 

PR.  GT.  .  551  ■  :• 

GO 

TO 

IF  CL. EG.  I 

.>  E I G-QPD I FF 

c  I 

IF  CL.  EG.  I 

.'  .  MnLL  =CPD  I 

FF  i  1." 

IF  CORD  IFF 

>:  I  >  .  GT  .  7  I G  • 

E I  G=DP  DI  FF  •.  I 

IF  CQPD IFF 

■I' . L  T . SMALL 

:•  SMALL  =  OPDIFF  C 

n 

SUM=3UM+0 

FI' IFF  .  :  • 

SUM 30=. SUM 

SO c-OPD  I  FF  •:  I  > 

♦  OP D IFF  •:  I 

i  ,i 
■~4  'Z,  'U 


GO  TD  4000 

3999  NDIV=NDIV  -  1 

4000  CONTINUE 

ShDPD-  >: SUM  ND I  V>  ♦  ■:  SUM.- ND I  V> 
ft  S  0  P  D =  S  U  M  S  O'  N  I*  1 V 
TI.IDP  I  =3 .  14]  ‘393654 

TMOP I  $0=1  NCR  I  ♦  TMGP  I 
MR  V  E  P  S  9  =  1.19  V  E  P  <• ! .  i  A  v  E  P 
R  M  S 1. 1 = h  S  □  D  -  ;  h  G  P  D 

STPEHL  =  1  .  -  '■  T U u P I S 0  l * l A V E P S 0  •>  *RMSO 

IF <STREHL.LT.  0.  •'  STPEHL  =  EXP  >:.  -  ■:  TWOP I S  0  -■  l.U-iVLT  PS  O’)  ♦RMS  GO  . 

D I  FF= BIG-SMALL 

♦  PRINTING  OUT  THE  FINAL  VALUES  AT  THE  RECEIVER  PLANE 

WRITE  CA.  504)  POUEP .  HAVE ■>  RANGE 

5 i_i 4  FORMAT  .  ♦  1  ♦ .  4  - ■, «  ♦ - — - —  —  — - — - 

$ - ♦  ./-••<  GX.  ♦  INDEX  VALUES  APE  FROM  PR  DP  MB  CODE  WHERE !  ♦  * 

iliz/c* UX .  ♦THE  BEAM  F’O'.lEF'  .  ♦  *  ol  U.  ’A  SX  ♦I.IhTTS^-*'X 
$.  c  U-9  ♦THE  l.iAVt LENGTH  HA..-.  ^ .  blc  .  O'  cX*  ♦  M I  lF'CN  S  ♦••**  •••* 

•I .  ci  1.1  X  '♦THE  P*  h  N  i E  W  h  _  ♦ .  A  X .  in  1  c  •  G .  cl  X .  ♦  M  E  T  E  P* .  ♦  / 

UPITE  5  03') 

503  FORMAT <3 OX- ♦ - ALL  MEASUREMENTS  OF  DISTANCE  IN  THIS  PROGRAM 

$  ARE  IN  .■•••■' METERS-"-'' - 

IF  iHPPINT.EA. O  '  GO  TO  8100 
IF CNPRIMT. EC. -1 >  GO  TO  8050 
DO  3000  1=1. NR AY 
WRITE CA. 505 1  I 

IF  >PAY '  I . 7> . LT. 0.  '•  GO  TO  7000 

48 


505  FORMAT  .  4X-  ♦  FOP  FRY*-' 13 

WRITE"  1  6  -  5  06')  NO  ••  I'>  »  ALPMO  'T-'  *  PAY  •'  1  «  1 ')  »  PAY  C I  -  £  > 

Mr'  I T  E  1  6 «  507  >  YD  <  I  ■•  -  hlFYD  <  I  >  -  PAY  <  I  *  3")  *  POV  •"  I  •  4 ) 

I  IP  1 7 1  •  6  •  5  03  >  F  RY  ■  1  •  5‘> 

503  FCr  • v’ :  •  3’ ■ .  o  XO  =*-613.  1£«4X«  *•  ALPHA  XO  =*»613.  12*4X 
3-  *  i-  ■  e  «  6 1 3 .  l  3 -  6:  '■  -  <-A*_ P HR  XF  =  *  -  5 1 3 .  1  cV  * 

50?  rri:  t;.-  -  yO  =*- 61 3.  15«4!Y*  ALPHA  YD  =*■*315:.  13?  4X 

1 «  »  i  F  ~  * «  6 1  3 .  12-  3: :  •  ♦ALPHA  YF  =  * «  5  1 3 .  3  £■■■'  > 

503  "I  ?.3' ■■'.«  *ALPHA  Zr  -*  -  51 3.  1 £•■' ') 

iiiR  I  TE  '  o -  5 O't* •'  PA  f  ■.  I  >  r  > 

505  FORMAT  <67X*  ♦  THE  FINAL  INDEX  DF  REFRACTION  IS  *  -  6 1 3 .  1  cy"> 

l, IP  I  T  E  • 6  -  5 1  0  ’  DP D I FF  •'  I 

510  FDPHHT  <6?X-  *  THE  TOTAL  OPTICAL  PATH  D I FFEPENCE  IS  * -  613 .  1£  •••'•■"'  •' 

GO  TO  3000 

7  0 0 0  IF  <'NPLAN  •:  I  .  E 0 .  1  GO  TO  7  1  0 0 
NP MINUS  =  NPLAN  ■  I  :■  -  1 

l, IP  I  TE  6-511  MO  I  .•  .  ALPXO  •"  I  )  -  NPM I NU  S  -  NPLAN  <  I "> 

511  FORMAT  <3X  -  «►  Xu  =  *■- 613.  1  3- 4X- ♦  ALPHA  XD  =*-613.  12-  4X 

3-  ♦THE  RAY  LEFT  THE  BEAM  OF  1 11  BETWEEN  PLANES  ♦« 13* ♦  AND*-  13*  *  •  <t  y 
MR  I  TE  >  6  -513  >  YO  1  I  •’  -  ALF'YD  I  > 

512  FORMAT <8X-  ♦  YO  =<>- 613.  12- 4X* *  ALPHA  YD  =  <  «  61  3 .  1 
GO  TO  3000 
7100  CONTINUE 

WRITE  <6-  5 1 3>  XO  <I)  -  ALF’XO  '-I  •'  -  NPLAN  <  I  > 

5 13  F D R M A T  3 X,  *  X 0  =  •«  61 3 .  1  2  -  4 X «  «•  A L P H A  X 0  =  «  -  6 1 3 .  1  £  *  4 X . 

3  *  T  H  E  PAY  LEFT  THE  BEAM  GRID  BETWEEN  THE  ThRhET  F'LhNE  ANB  PLANE*-  I 

3-  ♦  .  *■••) 

I..IP  I  TE  <6 -514'  YD  <  I  ')  -  ALPYD  <  I  ■' 

514  FORMAT  C3X-  ♦  YO  =«  -  G 1 3.  1  £- 4N- *  ALPHA  YD  613.  l£-'-'/> 

3000  CONTINUE 

8050  WRITE  < 6-515')  BIFF-  ASOPB-  SAOPD- STPEHL 

515  FORMAT  <25X* ♦THE  MAXIMUM  DEVIATION  ACROSS  THE  OUTPUT  PLANE  IS  ♦- 
-I'  G 13.  1  £  /  y  -  £  5  X  -  <•  O  R’  D  S  i.- 1 B  A  P  =  ** .  w  1 8 .  1  l  -  6  X  -  *  G  P  P  B  h  R ...  i.<  —  o  -  b  1  ,  1  c.  .-'  c  5 1  - 


3*THE  STPEHL  PA 
IF  f'MASti .  EC'.  £0> 
516  FORMAT <25X* ♦ — 

TIO  FOR  THIS 
WRITE <6- 516) 
- NOTE: 

PUN  IS* 

THIS 

-  613. 3  •'"  •••"  ) 

STPEHL  RAT  ID 

IS  MASK E B 

TD 

0 .  £ 

%  MET EPS. - 

IF  CM ASK . EC. 30 > 
IF  ■'MASK.  EC.  1  0> 
517  FORMAT  <£5X-  ♦— 

- <>  /  /  j 

1'IR'ITE  '.6*  517') 
WRITE <6- 513) 
- NOTE: 

THIS 

STPEHL  PATIO 

IS  MASKED 

TO 

0 .  3 

3  METER'S. - 

IF  (MASK.  EC.  4  0 :■ 
5££  FORMAT  <'£5X-  ♦ — 

- <>  X  /  j 

UPITE»:6^3£> 
- MOTE : 

THIS 

STPEHL  RATIO 

IS  masked 

TD 

0.4 

3  METER'S. - 

IF  CM  ASK.  EO.  5  O') 
523  FORMAT  <£5X-  «- — 

FIR ITE  <6-  5£3') 
- NOTE: 

7  H I S 

STPEHL  PATIO 

IS  MASKED 

TO 

0 . 5 

3  METERS. - 

IF  < MASK. EC. 55) 
524  FORMAT '£5X. ♦ — 

- 

MP I  TE  OSt524> 

- NOTE: 

THIS 

STPEHL  PATIO 

IS  MASKED 

TD 

0 . 55 

3  METERS. - 
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5 IS  FORMAT  i£5X»  ♦ - NOTE:  THIS  STPEHL  PATIO  IS  NR' TED  TO 

S-HETEPS  . - *'-') 

N  N  F  L  A  G  =  M  F  L  A  G  ~  1 

IF  >. HELPS.  RT.  1  :>  I.IP I TE  •  G <  5c 1  1  NNFLRG 
5£1  FGPMA’i  •Sr.:-:.  *«iss:j:;:;:THIS  TP  PS  E  MAS  DONE  WITH  THE  F  1 P  S  T*  <  3  S  •  c!  ; 
i.  PEL  r  *  .  ^  I  P  P‘E  D  .  •» ;;  *  ■  ■'  •' 

I F  f  i  HP'  I RT  .  EO.  1‘ .  CP.  •  NC'R  I  NT  .  EC1 .  - 1  ')  GO  TO  9000 
100  I  i*:'  1  T  ~  «  ’•SIP1  Pm.iEf  -  mihVEi  ' ?  F  httbi  i  HB_  ?  Mli  I V 
51'=*  Fn:'  ,-  '  l  BE  10.4*  £14) 

DU  SEOO  I  =  1  •  NR  AY 

♦  IF  PAY  *:  I  •  7  >  IS  LESS  THAN  0.  .  THhT  PAY  IS  BYPASSED. 

IF  (PRY  I  ?  ?')  .  LT  .  0.  :*  GO  TO  Sc  0  0 
□PBIF-Or-TiIFF  >1  >  .-•'I, l AVER 

MR  1  T r.  1  'P  <  5c  U  >  PAY  ( 1 5  1 )  i  PAY  ( I »  S'  •*  «  OF' D  I F 
5c  0  FORMAT  •  5S  <  SF  1  0 .  4  :■ 

S00  CONTINUE 
000  STOP 
ENB 


SUBROUTINES  BEGIN  HERE 


SUBROUTINE  XYNEM  CKDLB*  ALPHhSH  YOLlH  ALPHR'tU  Zf-'OLB'  ZF'NEXT  ,  ALt^R 
SYNE  W) 

<>  CALCULATE  NEW  X  Y  lOGRBINhTEc- 

XNEW  =  XOLD  +  CZRNEXT  -  ZPOLB)  •  ALPHAS-- ALPHAS') 

YNEW  =  YOLB  +  t'ZPNEXT  -  ZF'OLD)  *  (  ALPHAYVALPHAZ') 

RETURN 

ENB 


SUBROUTINE  AVAN  (RNZEPO*  hNONE  «  Art  FAR') 

♦  CALCULATE  AVERAGE  INDEX  BETWEEN  MIDFLAMES 

AN  BAR  =  <  ANZEPD  +  ANDNt.)  --’c 

RETURN 

ENB 


sue  pout  i  nr  opl  -ope.  xold«  xhew.  vole,  yhem.zdld.  znei.u»'  •  •  hp> 

*  CALCULATE  OPTICAL  PATH  LEMETH  BETWEEN  MIHELhHE. 
eel::  -  XHEI.I  -  XCLB 
EELV  =  MEM  -  VOL  D 
I'ELI  -  IHF.M  -  IDLE 

p  s  :opt  •  belx*eelx  +  belt ♦eel v  +  delz*delz:> 

OPE  —  OPE  +  HtiFHf'  *P 

PETUPM 

END 


SUBPOUT  1  HE  IP  POP  •  r:CLE«  ISUl.H  EELZMI'H  ZPNU> 

♦  FI  ME  HIli1  LOCATIONS  Or  PH  Ft  5  E  SCPEEHS  AMD  THE  IP 

♦  PEI  PEC T I VE  MIEPLENEI 
2PNU  "  HOLE  +  ■  EELItfi.L- 
2SNU  =  HOLE  +  EELZHi.l 
PETUPM 

EME 


SUBPOUT  I  HE  ANCLE • TZ I N. THETA) 

♦  F I  HE  ANCLE  OP  PHY  TO  SUPFhi.E  HOF’MHL 

♦  TZIM  =  EIPECTION  C 01  I  HE  OF  PHY  WITH  Z  AXIS 
THET  fl”HCOS (TZ I N) 

PETUPM 

END 


SUBPOUT  I  HE  PFIJM I  TV  (' TX  I M  <  TV  I H  <  T  Z I M »  HM I H  -  HNPF  <  THETAP) 

♦  FI HE  UNIT  VECTOR  OF  THE  PEFPAC7EB  PHY 

♦  txih.  TYirn  *,  tzim  ape  bipectioh  cosines  of  imcieemt  phy 

♦  DM  IMPI.IT  HUE  THE  EIPECTION  COSINES  OF  THE  PEFFACTED  PHY 

♦  OH  OUTPUT 

♦  AMIN  (HNPF)  IS  THE  INDEX  OF  THE  IMCIEEMT  (P  ESP  ACTED')  PAY 

♦  THETAP  =  PEFPftCT ION  ANCLE 

THE T HP = AS  I M  <  ■  h r 1 1  N-HNr  F)  I H  ■  T  HETRP)  > 

T X I N -  (  AH  I H , •  R  H  F  F ')  ♦  T X I H 
TYIM=  (RN  IN.- HNPF)  «-TYIH 
TZ IH=C OS  (THETHP) 

PET UP H 

END  51 


♦  ♦  «  <  <■<  *  <  *  '  <•  <•  «■  *■ 


s  utf  n  i  me  ru  r  z  ■  alpha's  .  alf-hay  ,  alphaz> 

♦  FIND  ALPHA  Z  PGP  INPUT  VALUE i  OF  ALPHA  X  ALPHA  V 
alphas  -  i opt  •.  i .  -  alrhaxoalphax  -  alphay ♦alphay » 

crjlj:  :i 
END 


S U I* fr‘ OUT  1 N tf.  h N 1 1 E :  •  '•  I* E  L X * i’ « X P ft Y  «  r  P ft Y  ♦  ML  *  E M D E  F « H I  M L E X 
Ii  I M E N s  ION  X  <  N  L  >  *  Y  •  N L>  >  A I NDE X  <NL»  NL > 


<*  C-  O  C>  <•  V  (►  <>  < 


V-  IPANIJM) 


♦  PROPMD  SCREEN  S  Ac  E  CALCULATEB  1.1 1 TH  RESPECT  TO  THE  XMITTEP 

♦  PLANE  AND  AP  PAY  VALUE  -1«1>  IS  IN  THE  UPPER  LEFT-HAND  CDF  HEP 

♦  OF  THE  GRIB  WHEN  LOOPING  FROM  THE  XMITTEP.  THEREFORE.  MY 

♦  INTERPOLATION  IS  FROM  THE  UPPER  RIGHT-HAND  CORNER. 

♦  XRAY  =  X  GRIB  COORDINATE  OF  THE  RAY 

♦  YRAY  =  Y  GRID  CODRDINhTE  Dr  THE  FAY 

♦  NL  =  ARRAY  SIZE  OF  THE  RPOPND  PHASE  SCREEN 

♦  DELXY  =  GRID  SPACING  OF  THE  PR OP MB  PHASE  SCREEN 

♦  IRANIJM  =  NUMBER  OF  THE  RAY  HI  THIN  ITS  BUNDLE 

♦  I  XL  IXH  =  X IXL>  >  XRA'y  >  Y  1  I XH > 

♦  IYL  C  IYH  =  Y  <.  I  i'L  ■'  YRAY  Y  <  I YH1 

♦  ENDEXF  =  FINAL  INTERPOLATED  VALUE  OF  THE  INDEX 


IF  U RANUM . GT . 1 >  GO  TO  1 0 0 

♦  ONLY  WANT  TO  SET  UP  X  C  Y  MATRICES  FOR  FIRST  RAY  OF  THE  BUNDLE 


♦  SET  UR  COORDINATE  GRIDS  IN  X  %.  Y  DIRECTIONS 

♦  WITH  CENTER  OF  GRID  •  0,0:.  BTWH  X  <'NL£>  &  X'ULE  +  l) 


DEL  X  Y  S = D  E  L  X  Y / 2 
NL2  =  NL-'E 
DO  99  1=1 -NL 
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150  XP  =  XRRY  -  X  < I XL> 

VO  ~  YPRY  —  Y  '■  1  i  L  ■' 

ENDEXL  -  H INDEX < 1YL«  IXL  >  + 
*DFLXY>  ♦XP 

ENDEXH  =  Hi  NDEX  •'  IYH«  I  XL>  + 

•i-I'iE-il:  :  v  :■  «xp 

ENDEXF  =  ENDEXL  +  ■ ENDEXH 

RETURN 

£00  C  UNT  I HUE 

*  PRY  COORDINATES  hPE  OUTS 
EHDEXE  =  -900. 

RETURN 

END 


:  >:  n  I NDEX  <  I YL  »  I XH  > -ft I NDEX  <  I YL «  I  XL  ">  > 
..Hi  riDEX  1  I YH .  I  XH  •  -H 1  HDEX  <  I  YH«  I  XL  >  > 
EHI'EXL  :>  /DELXY>  ♦YO 

iE  GPIDBOUNDS — SET  INDEX  NEGATIVE 
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PPD3P  hM  I.  AT  ft  I  riPI.IT «  OUTPUT  *  TAPE88  *  T  APE8  1  "•' 

D I HEM1  1  CM  MPLhME  <  1 1  '•  «  DELE  <  1  l'1  *  ft  •  64 »  64 >  «  P  •  64 « 64')  ,  C  < 64 *  64 
D 1  ME M  7  I D U  D  *  6. 4 »  6-4 )  «  E  1  6 4 »  64  ■  ,  F  •.  6-4 «  6-4  )  «  Ii E  L  7 ' Y  <  1  1  1  «  D  E  L  T  ft _ 
Ii  I  HE i  ■  7  I  I’M  F  P.4 »  P'4  >  ,  H  •  64  ♦  64  >  «  P  1  P.4  •  64  >  ,  :  «.  6-4  •  64  >  «  T  <  64 ,  64  > 

F E ft7'-  1  c  I* 1  POMEF,  l.lftVE *  F'MMitc.  ,  i~'ETftM,  MDZ 

CFftl  '88  •  DELXY  <  1  ■-  ,  MLM 

F‘E ("7 •  ■  -  '  Mi-'LftMt;  •  1.)  ,  D  — .  L  4  <  1  ) 

C  P  I'L  ’  “•  c'  ♦  -! .  1  4  1  5 lilH'v'E 

DC  9  :«  I  =  1  <  MLM 
DO  96  .i=i,  MLM 
•F* Pi  T  •  I ,  J  1  =  1  ,  ij  +  PEI  hM 
99  CDMTIMUE 

PERU  1  f  c  >  DEL X Y  <  8  -' 

HD  5  1=1 5 MLM 

5  PE  AH  <'88')  <  Z j  ,  j.<  ,  j=  l  ,  ml  M) 

PE r:D  <88>  MPLhME  <£'>  >  DELZ  <£> 

P  E  ft  D  <  8  8  )  DEL  X  V  <  8  > 

DO  6  1=1? MLM 

6  PE  ftli  <  88  >  <  P  <  I ,  J  >  ,  J=  1 ,  MLM  ‘ 

PEftD  <'££'>  MPLftME  <3>  »  DELZ  <'3> 

P E  ft D  <'££■>  I'EL ■■  V  < 4  ) 

IiQ  7  1  =  1  J  MLM 

7  P  E  ft  Ii  <  8  £  ">  <■  H  <  I  >  J  7'  ?  J  =  1 ,  M  L  H  "> 

PEftD  <887*  MPLftME  C4>  ,  DELZ  <4*> 

PEftD  <  33>  DEL  XV  <  5') 

DO  8  1  =  1, MLM 

8  PEA  D  cl  cl  ">  <  8  <'  I  *  J  >  ,  J=1  ,  M  L  M  '> 

PEftD  <££>  MPLftME  -  5>  ,  DELZ  <57' 

PE  AD  <88<DELXY<  67. 

DO  9  1=1, MLM 

9  P  E  ft  D  <  8  8  7>  <  F  ('  I ,  J  7*  , .  J  =  1 ,  M  L  M ') 

PEftD  <887>  MPLftME  <6'»  »  DELZ  <'6> 

PEftD  <  88 1 DELXY <  <  > 

DO  10  1=1, MLM 

I  0  PEftD  <  88>  c'E  <  I  ,  J 7'  ,  J=  1 ,  MLM') 

READ  <887'  MPLftME  < 7),  DEL 4  <7) 

PEftD  <88 ')  D  E  L  X  V  <  8  '■' 

DO  11  1=1 5 MLM 

I I  PEftD  <887'  <D  <  1 ,  J7'  ,  J=  1  ,  MLM 7* 

PEftD  <  8  8  :•  M  P  L  ft  M  E  <  8  7>  ,  DELZ  <  3  > 

READ  <887'  DELXY  <97* 

DO  18  1=1, MLM 

18  PEftD  <  88  >  <C  <  I ,  J  1  ,  -l=  1 ,  MLM) 

PEftD  <887'  MPLftME  <9>  «  DELZ  <9> 

P  E  ft  D  <  8  8  7*  DEL  X  Y  <  1  U  > 

DO  13  1=1, MLM 

13  PEftD  <  88>  <  E  <  I ,  -  J  ■  ,  .J=l  ,  MLM) 

PEftD  <88>  MPLftME  <1  O')  ,  DELZ  (10') 

PEftD < 88>  DELXY  <11 ) 

DO  14  1=1, MLM 
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F  EhU  •  ££'>  •'  ft  1  I  *  J>  •  J=  1  *  NLM') 

PF.h l'i •  f  c  ■  HPLftNE  ■'  11*-  DELZ  ■  1 1> 
F'ftf'-  o  -  F'ftNbE 
do  1 1.  i  =  i  •  1 1 

DEL'!  r* Z  1  1  >  --  Phi (i-  —  DELZ  <  lc'  —  I  1 

Riv-c-  =  delz  •  1  e-  I  > 

I F  ■  I  .  E 0 .11)  ED  T □  IE 

Cu'.'i  I Ml.iE 

CONTINUE 

DO  E'E'  1  =  1  .NLM 

ft  *■  1  •  1'j  =  0 . 

E  <  I »  1 )  =  0. 

C  '1-1  =  0 . 

D  <  I »  1  .)  =  Ci . 

Ea»ii=  o. 

F  C I »  1 )  =  0. 

0  ■ I  j  1  >  =  0 . 

H  C  I  -  1  '>  =  0 . 

P  <1,1 > =  0 . 

S  •:  I »  1  :■  =  0 . 


DO  El  J 

=  1  *  NLM 

ft  <  I «  J >  = 

ft  •:  I  *  J> 

-Cl  LIES 

+ 

TCI* 

E  <  I » -.1  >  = 

E  ■:  I  *  J) 

••'Cf:  del: 

+ 

T  CI  , 

CCD  J  :•  = 

C  C  I  *  J  - 

C  K  l.i  E  X 

-t- 

TCI, 

I'  <  I  *  J  )  = 

D  <  1  *  J> 

•'CKDEX 

+■ 

T  Cl, 

E  CI  * . J )  = 

ECDJ) 

CKDEX 

+■ 

TCI* 

F  C I ,  ..l  )  = 

F  ■.  I  *  J  :• 

••'Ci  del: 

+ 

T  •:  I  * 

G  •  I  ,  j  f  - 

G  <  I  *  J ) 

••'CKDEX 

+ 

T  'I  * 

H  C I  >  J  :■  = 

H  -  I ,  J  ' 

••'L-K  DEL’ 

+ 

TCI, 

Pci*  J )  = 

PCI*  J  :• 

•ck  del: 

+ 

TCI* 

S  < I  9  J)  = 

S  C  I  *  -J ) 

••'CKDEX 

+• 

TCI, 

CONTINUE 

WRITE  (£1.)  F'OWiP*  WftVE >  F'ftNi-E*  F'ETftM 
♦  FIRST  F'LftHE  IS  ThPGET  F'LRNE 
WRITE':' El  5  NDZ*  NLM 
WRITE  CEl  )  DELTftZ  C  1  ') 

WP ITE  *:  E1  N P L ft H E  <  1>  *  DEL XV  <  1 1  ’> 

DO  31  1=1 *NLM 


WRITE  Cel 

*:  ft  <  I  *  J 

*  J= 1  *  NLM  ) 

WR I TE  C£ 1 

DELTftZ 

C  £  ) 

WRITE  Cel 
DO  3£  1= 

NPLfiNE 

1  *  NLM 

'£.'*  DEL  XV 

<  1  0  > 

WRITE 'El 

•E • I  *  J 

*  J= 1  *  NLM) 

MR I TE  C£ 1 

DELTftZ 

»■’  ‘j  > 

MR  I TE  C£ 1 

NPLftNE 

C3) *  DEL XV 

DO  33  1=1* ULM 

WR I TE  C£  1  ')  •:  C  I  *  J  >  *  J=  1  *  NLM') 

WRITE  '-El'  DELTftZ  •' 4> 
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3*3 


Of 


WR  I  TE  •:  3  D  NRLRNE  (4  >  ,  DELXY  <8> 

I'D  34  1 1  «  NLM 

WR  I T E  •  3  1  •  •  B  >  I «  J  ■  « .J=  1 ,  HL.H> 

1.1 RDE  f£i  DELThZ  <5:> 

Wi;  1 1  “  •'  3  1  >  NRLhH£  <5  :•  ,  DELXY  *.?') 

Dll  1  —  1  f  HLM 

WR  I  ’  E  ■  3  1  ■  >  £  >  I ,  J  :<  ,  .J=  1 ,  HLt  V> 

Ur-1  ,  j;  •  c.l  >  DELThZ 

WR  L  T L"  1_'  1  ■'  HPLhME  1  E  '  ,  DEL XV  *:.6> 

DO  33  I  —  1 «  HLM 

WR  I  111  'c  1  ■'  ■'  r  ■  I  «  J  >  <  J=  1  i  HLM  •' 


WRITE 

WRITE 


■  3 1  )  BELT  R2  <  7  :< 

3 1  :•  URL  HUE  *.  ?':•  -  DEL XV  (5') 


DO  37  1=1, HLM 


WRITE 
WRITE  •: 
WRITE  •: 


C  1  )  '•  H  *.  I ,  .J  .*  ,  .J=  1  ,  HLM  1 

3 1 >  DELThZ •: 3 ' 

3D  HPLRHE  «:  *')  ->  DELXY  <4'» 


DO  38  1=1, HLM 


WRITE (31 
WRITE  (31 
WR I TE  (3 1 


J=  1  «  HLM':* 

, DELXY (3) 


40 


41 


i.  H  *  I ,  J  > 

':*  DELThZ  (3 
:>  MPLhME  O' 

DO  39  1=1, HLM 

WR  I  TE  (3 1  :•  (P  ■  I  *  J  >  »  J=  1 »  HLM':* 

WRITE  (31'>  DELTRZ  *1  0> 

WR  I  TE  (3 1  >  MPLFiME  '  1  O’)  ,  DELXY  (£> 
DO  40  1=1 , HLM 
WR  I  TE  (3 1  *.  3  ( I ,  J>  ,  J=1 ,  HLM  > 

WRITE  (3D  DELTRZ  dl> 

WR  I  TE  (3 1  >  HRLhHE  ( 1  1  *  ,  DEL XV  <  1  '• 
DO  41  1=1, HLM 
WR  I  TE  ( 3 1  >  (T  ( I ,  S>  >  J=  1 ,  HLM) 
PRINT*, ”  ” 

PR  I  NT  * , POWER , WhVE , RfiNGE 
PRINT*,"  " 

PRINT  *,  (R  (33, J> , J=1 , 1 0> 

ENDFILE  31 

STOP 

END 


ViXA 
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